# yI0-xI0.R
# Philips, Andrew Q. "How to Avoid Incorrect Inference (While Gaining Correct Ones) in Dynamic Models". Forthcoming at Political Science Research and Methods. 
# andrew.philips@colorado.edu
# 2/9/21
#
# Figures produced in this R script:
#   --Main document, Figure 1: "yi0-xi0-SR.pdf"
#   --Main document, Figure 2: "yi0-xi0-LR.pdf"
#   --SI, Figure 1: "yi0-xi0-SR-evar5.pdf"
#   --SI, Figure 2: "yi0-xi0-LR-evar5.pdf"
#   --SI, Figure 3: "yi0-xi0-SR-MSE.pdf"
#   --SI, Figure 5: "yi0-xi0-LR-MSE.pdf"
#   --SI, Figure 4: "yi0-xi0-SR-MSE-evar5.pdf"
#   --SI, Figure 6: "yi0-xi0-LR-MSE-evar5.pdf"
# -----------------------------------

rm(list = ls())
setwd("~/Dropbox/My_Folder/Papers-Projects/Wlezien-Enns-PSRM/Final-Feb2021/scripts/")
# -----------------------------------
library(dynamac)
library(nlWaldTest)
library(ggplot2)
library(grid)
library(gridExtra)
library(dplyr)
library(directlabels)

set.seed(098765)
# -------------------------------------------------------------
# EXPERIMENT 1: Y ~ I(0), X ~ I(0)
sims = 2000 # no. sims
ar.param <- seq(0,0.95, by = 0.05) # autoregression in Y, X
N <- c(50, 250, 1000) # number of obs
e.var <- c(1, 5) # error variance in Y

container.ardl <- matrix(NA, nrow = sims*length(ar.param)*length(ar.param)*length(N)*length(e.var), ncol = 11) # null matrix to hold output
# cols = sims, ar.param.X, ar.param.Y, N, Beta.x, UL.beta.x, LL.beta.x, LR.x, LR.x.UL, LR.x.LL, e.var
container.ecm <- matrix(NA, nrow = sims*length(ar.param)*length(ar.param)*length(N)*length(e.var), ncol = 11) # null matrix to hold output
# cols = sims, ar.param.X, ar.param.Y, N, Beta.x, UL.beta.x, LL.beta.x, LR.x, LR.x.UL, LR.x.LL, e.var
container.ldv <- matrix(NA, nrow = sims*length(ar.param)*length(ar.param)*length(N)*length(e.var), ncol = 11) # null matrix to hold output
# cols = sims, ar.param.X, ar.param.Y, N, Beta.x, UL.beta.x, LL.beta.x, e.var
container.static <- matrix(NA, nrow = sims*length(ar.param)*length(ar.param)*length(N)*length(e.var), ncol = 8) # null matrix to hold output
# cols = sims, ar.param.X, ar.param.Y, N, Beta.x, UL.beta.x, LL.beta.x, e.var
# -------------------------------------------------------------

row <- 1
prog <- txtProgressBar(min = 0, max = sims, style = 3)
for (s in 1:sims) {
  for (ar.y in ar.param) {
    for (ar.x in ar.param) {
      for (n in N) {
        for (ev in e.var) {
          setTxtProgressBar(prog, value = s)
          # create Y
          if (ar.y == 0) {
            y <- rnorm(n + 100, sd = sqrt(ev))
          } else {
            y <- arima.sim(n = (n + 100), list(ar = c(ar.y)), innov = rnorm(n = (n + 100), sd = sqrt(ev)))
          }
          # create X
          if (ar.x == 0) {
            x <- rnorm(n + 100)
          } else {
            x <- arima.sim(n = (n + 100), list(ar = c(ar.x)))
          }
          # ARDL model:
          res.ardl <- lm(y[101:length(y)] ~ lshift(y, 1)[101:length(y)] + x[101:length(x)] + lshift(x, 1)[101:length(x)])
          # ECM model:
          res.ecm <- lm(dshift(y)[101:length(y)] ~ lshift(y, 1)[101:length(y)] + dshift(x)[101:length(x)] + lshift(x, 1)[101:length(x)])
          # LDV model
          res.ldv <- lm(y[101:length(y)] ~ lshift(y, 1)[101:length(y)] + x[101:length(x)])
          # static model
          res.static <- lm(y[101:length(y)] ~ x[101:length(x)])
      
          # grab QOI (FORALL): ----
          container.ardl[row, 1] <- container.ecm[row, 1] <- container.ldv[row, 1] <- container.static[row, 1] <- s
          container.ardl[row, 2] <- container.ecm[row, 2] <- container.ldv[row, 2] <- container.static[row, 2] <- ar.x
          container.ardl[row, 3] <- container.ecm[row, 3] <- container.ldv[row, 3] <- container.static[row, 3] <- ar.y
          container.ardl[row, 4] <- container.ecm[row, 4] <- container.ldv[row, 4] <- container.static[row, 4] <- n
          container.ardl[row, 11] <- container.ecm[row, 11] <- container.ldv[row, 11] <- container.static[row, 8] <- ev
        
          # QOI (ARDL)
          container.ardl[row, 5] <- coef(res.ardl)[3] # beta.x
          container.ardl[row, 6] <- confint.default(res.ardl, parm = c(3), level = 0.95)[2] # UL
          container.ardl[row, 7] <- confint.default(res.ardl, parm = c(3), level = 0.95)[1] # LL
          # LR effects (ARDL)
          lrm <- nlConfint(res.ardl, "(a[3] + a[4])/(1-a[2])", level = 0.95)
          container.ardl[row, 8] <- lrm[1] # LRM
          container.ardl[row, 9] <- lrm[3] # LRM UL
          container.ardl[row, 10] <- lrm[2] # LRM LL
      
          # QOI (ECM)
          container.ecm[row, 5] <- coef(res.ecm)[3] # beta.d.x
          container.ecm[row, 6] <- confint.default(res.ecm, parm = c(3), level = 0.95)[2] # UL
          container.ecm[row, 7] <- confint.default(res.ecm, parm = c(3), level = 0.95)[1] # LL
          # LR effects (ECM)
          lrm <- nlConfint(res.ecm, "(a[4])/(-a[2])", level = 0.95)
          container.ecm[row, 8] <- lrm[1] # LRM
          container.ecm[row, 9] <- lrm[3] # LRM UL
          container.ecm[row, 10] <- lrm[2] # LRM LL
      
          # QOI (LDV)
          container.ldv[row, 5] <- coef(res.ldv)[3]
          container.ldv[row, 6] <- confint.default(res.ldv, parm = c(3), level = 0.95)[2] # UL
          container.ldv[row, 7] <- confint.default(res.ldv, parm = c(3), level = 0.95)[1] # LL
          # LR effects (LDV)
          lrm <- nlConfint(res.ldv, "(a[3])/(1-a[2])", level = 0.95)
          container.ldv[row, 8] <- lrm[1] # LRM
          container.ldv[row, 9] <- lrm[3] # LRM UL
          container.ldv[row, 10] <- lrm[2] # LRM LL
      
          # QOI (static)
          container.static[row, 5] <- coef(res.static)[2]
          container.static[row, 6] <- confint.default(res.static, parm = c(2), level = 0.95)[2] # UL
          container.static[row, 7] <- confint.default(res.static, parm = c(2), level = 0.95)[1] # LL
      
          row <- row + 1
        }
      }
    }
  }
}
head(container.ardl)
head(container.ecm)
head(container.ldv)
head(container.static)
#save.image("scenario-yi0-xi0-data.RData")

# Create rejection rates and label:
container.ardl <- as.data.frame(container.ardl)
# cols = sim, ar.param.X, ar.param.Y, N, Beta.x, UL.beta.x, LL.beta.x, LR.x, LR.x.UL, LR.x.LL, e.var
colnames(container.ardl) <- c("sims", "ar.param.X", "ar.param.Y", "Time", "Beta.x", "Beta.x.ul", "Beta.x.ll", "LR.x", "LR.x.UL", "LR.x.LL", "e.var")
container.ecm <- as.data.frame(container.ecm)
# cols = sim, ar.param.X, ar.param.Y, N, Beta.x, UL.beta.x, LL.beta.x, LR.x, LR.x.UL, LR.x.LL, e.var
colnames(container.ecm) <- c("sims", "ar.param.X", "ar.param.Y", "Time", "Beta.x", "Beta.x.ul", "Beta.x.ll", "LR.x", "LR.x.UL", "LR.x.LL", "e.var")
container.ldv <- as.data.frame(container.ldv)
# cols = sim, ar.param.X, ar.param.Y, N, Beta.x, UL.beta.x, LL.beta.x, LR.x, LR.x.UL, LR.x.LL, e.var
colnames(container.ldv) <- c("sims", "ar.param.X", "ar.param.Y", "Time", "Beta.x", "Beta.x.ul", "Beta.x.ll", "LR.x", "LR.x.UL", "LR.x.LL", "e.var")
container.static <- as.data.frame(container.static)
# cols = sim, ar.param.X, ar.param.Y, N, Beta.x, UL.beta.x, e.var
colnames(container.static) <- c("sims", "ar.param.X", "ar.param.Y", "Time", "Beta.x", "Beta.x.ul", "Beta.x.ll", "e.var")

# proportion false rejection
container.ardl$reject.beta.x <- ifelse(container.ardl$Beta.x.ll > 0 | container.ardl$Beta.x.ul < 0, 1, 0)
container.ecm$reject.beta.x <- ifelse(container.ecm$Beta.x.ll > 0 | container.ecm$Beta.x.ul < 0, 1, 0)
container.ldv$reject.beta.x <- ifelse(container.ldv$Beta.x.ll > 0 | container.ldv$Beta.x.ul < 0, 1, 0)
container.static$reject.beta.x <- ifelse(container.static$Beta.x.ll > 0 | container.static$Beta.x.ul < 0, 1, 0)

# coverage of LRM:
container.ardl$reject.lrm <- ifelse(container.ardl$LR.x.LL > 0 | container.ardl$LR.x.UL < 0, 1, 0)
container.ecm$reject.lrm <- ifelse(container.ecm$LR.x.LL > 0 | container.ecm$LR.x.UL < 0, 1, 0)
container.ldv$reject.lrm <- ifelse(container.ldv$LR.x.LL > 0 | container.ldv$LR.x.UL < 0, 1, 0)

# squared error:
container.ardl$sqerr.beta <- (container.ardl$Beta.x - 0)^2
container.ecm$sqerr.beta <- (container.ecm$Beta.x - 0)^2
container.ldv$sqerr.beta <- (container.ldv$Beta.x - 0)^2
container.static$sqerr.beta <- (container.static$Beta.x - 0)^2
container.ardl$sqerr.lrm <- (container.ardl$LR.x - 0)^2
container.ecm$sqerr.lrm <- (container.ecm$LR.x - 0)^2
container.ldv$sqerr.lrm <- (container.ldv$LR.x - 0)^2


# collapse everything down:
container.collapse.ardl <- container.ardl %>%
  group_by(ar.param.X, ar.param.Y, Time, e.var) %>%
  summarize(reject.beta.x = mean(reject.beta.x),
            reject.lrm = mean(reject.lrm),
            mse.beta.x = mean(sqerr.beta),
            mse.lrm = median(sqerr.lrm))

container.collapse.ecm <- container.ecm %>%
  group_by(ar.param.X, ar.param.Y, Time, e.var) %>%
  summarize(reject.beta.x = mean(reject.beta.x), 
            reject.lrm = mean(reject.lrm),
            mse.beta.x = mean(sqerr.beta),
            mse.lrm = median(sqerr.lrm))

container.collapse.ldv <- container.ldv %>%
  group_by(ar.param.X, ar.param.Y, Time, e.var) %>%
  summarize(reject.beta.x = mean(reject.beta.x), 
            reject.lrm = mean(reject.lrm),
            mse.beta.x = mean(sqerr.beta),
            mse.lrm = median(sqerr.lrm))

container.collapse.static <- container.static %>%
  group_by(ar.param.X, ar.param.Y, Time, e.var) %>%
  summarize(reject.beta.x = mean(reject.beta.x),
            mse.beta.x = mean(sqerr.beta))

# separate by T:
container.collapse.ardl.50 <- subset(container.collapse.ardl, Time == 50)
container.collapse.ardl.250 <- subset(container.collapse.ardl, Time == 250)
container.collapse.ardl.1000 <- subset(container.collapse.ardl, Time == 1000)
container.collapse.ecm.50 <- subset(container.collapse.ecm, Time == 50)
container.collapse.ecm.250 <- subset(container.collapse.ecm, Time == 250)
container.collapse.ecm.1000 <- subset(container.collapse.ecm, Time == 1000)
container.collapse.ldv.50 <- subset(container.collapse.ldv, Time == 50)
container.collapse.ldv.250 <- subset(container.collapse.ldv, Time == 250)
container.collapse.ldv.1000 <- subset(container.collapse.ldv, Time == 1000)
container.collapse.static.50 <- subset(container.collapse.static, Time == 50)
container.collapse.static.250 <- subset(container.collapse.static, Time == 250)
container.collapse.static.1000 <- subset(container.collapse.static, Time == 1000)

# plot: -------------------------------------------
prob <- c(.10, 0.25, 0.50) # breaks for contours

# Short-run:
p1 <- ggplot(data = subset(container.collapse.ardl.50, e.var == 1), aes(x = ar.param.X, y = ar.param.Y, z = reject.beta.x)) + theme_minimal() + geom_raster(aes(fill = reject.beta.x)) + ylab("Autoregression in Y") + xlab("Autoregression in X") + ggtitle("T=50, ARDL/ECM") + guides(fill=guide_legend(title="Proportion \n Rejected")) + scale_fill_continuous(limits=c(0, .7), breaks=seq(0, .7, by = 0.1), low = "white", high = "black") + geom_contour(aes(z = reject.beta.x), breaks = prob) + geom_dl(aes(label = ..level..), method="bottom.pieces", stat="contour", breaks = prob)

p2 <- ggplot(data = subset(container.collapse.ldv.50, e.var == 1), aes(x = ar.param.X, y = ar.param.Y, z = reject.beta.x)) + theme_minimal() + geom_raster(aes(fill = reject.beta.x)) + ylab("Autoregression in Y") + xlab("Autoregression in X") + ggtitle("T=50, LDV") + guides(fill=guide_legend(title="Proportion \n Rejected")) + scale_fill_continuous(limits=c(0, .7), breaks=seq(0, .7, by = 0.1), low = "white", high = "black") + geom_contour(aes(z = reject.beta.x), breaks = prob) + geom_dl(aes(label = ..level..), method="bottom.pieces", stat="contour", breaks = prob)

p3 <- ggplot(data = subset(container.collapse.static.50, e.var == 1), aes(x = ar.param.X, y = ar.param.Y, z = reject.beta.x)) + theme_minimal() + geom_raster(aes(fill = reject.beta.x)) + ylab("Autoregression in Y") + xlab("Autoregression in X") + ggtitle("T=50, Static") + guides(fill=guide_legend(title="Proportion \n Rejected")) + scale_fill_continuous(limits=c(0, .7), breaks=seq(0, .7, by = 0.1), low = "white", high = "black") + geom_contour(aes(z = reject.beta.x), breaks = prob) + geom_dl(aes(label = ..level..), method="bottom.pieces", stat="contour", breaks = prob)

p4 <- ggplot(data = subset(container.collapse.ardl.250, e.var == 1), aes(x = ar.param.X, y = ar.param.Y, z = reject.beta.x)) + theme_minimal() + geom_raster(aes(fill = reject.beta.x)) + ylab("Autoregression in Y") + xlab("Autoregression in X") + ggtitle("T=250, ARDL/ECM") + guides(fill=guide_legend(title="Proportion \n Rejected")) + scale_fill_continuous(limits=c(0, .7), breaks=seq(0, .7, by = 0.1), low = "white", high = "black") + geom_contour(aes(z = reject.beta.x), breaks = prob) + geom_dl(aes(label = ..level..), method="bottom.pieces", stat="contour", breaks = prob)

p5 <- ggplot(data = subset(container.collapse.ldv.250, e.var == 1), aes(x = ar.param.X, y = ar.param.Y, z = reject.beta.x)) + theme_minimal() + geom_raster(aes(fill = reject.beta.x)) + ylab("Autoregression in Y") + xlab("Autoregression in X") + ggtitle("T=250, LDV") + guides(fill=guide_legend(title="Proportion \n Rejected")) + scale_fill_continuous(limits=c(0, .7), breaks=seq(0, .7, by = 0.1), low = "white", high = "black") + geom_contour(aes(z = reject.beta.x), breaks = prob) + geom_dl(aes(label = ..level..), method="bottom.pieces", stat="contour", breaks = prob)

p6 <- ggplot(data = subset(container.collapse.static.250, e.var == 1), aes(x = ar.param.X, y = ar.param.Y, z = reject.beta.x)) + theme_minimal() + geom_raster(aes(fill = reject.beta.x)) + ylab("Autoregression in Y") + xlab("Autoregression in X") + ggtitle("T=250, Static") + guides(fill=guide_legend(title="Proportion \n Rejected")) + scale_fill_continuous(limits=c(0, .7), breaks=seq(0, .7, by = 0.1), low = "white", high = "black") + geom_contour(aes(z = reject.beta.x), breaks = prob) + geom_dl(aes(label = ..level..), method="bottom.pieces", stat="contour", breaks = prob)

p7 <- ggplot(data = subset(container.collapse.ardl.1000, e.var == 1), aes(x = ar.param.X, y = ar.param.Y, z = reject.beta.x)) + theme_minimal() + geom_raster(aes(fill = reject.beta.x)) + ylab("Autoregression in Y") + xlab("Autoregression in X") + ggtitle("T=1000, ARDL/ECM") + guides(fill=guide_legend(title="Proportion \n Rejected")) + scale_fill_continuous(limits=c(0, .7), breaks=seq(0, .7, by = 0.1), low = "white", high = "black") + geom_contour(aes(z = reject.beta.x), breaks = prob) + geom_dl(aes(label = ..level..), method="bottom.pieces", stat="contour", breaks = prob)

p8 <- ggplot(data = subset(container.collapse.ldv.1000, e.var == 1), aes(x = ar.param.X, y = ar.param.Y, z = reject.beta.x)) + theme_minimal() + geom_raster(aes(fill = reject.beta.x)) + ylab("Autoregression in Y") + xlab("Autoregression in X") + ggtitle("T=1000, LDV") + guides(fill=guide_legend(title="Proportion \n Rejected")) + scale_fill_continuous(limits=c(0, .7), breaks=seq(0, .7, by = 0.1), low = "white", high = "black") + geom_contour(aes(z = reject.beta.x), breaks = prob) + geom_dl(aes(label = ..level..), method="bottom.pieces", stat="contour", breaks = prob)

p9 <- ggplot(data = subset(container.collapse.static.1000, e.var == 1), aes(x = ar.param.X, y = ar.param.Y, z = reject.beta.x)) + theme_minimal() + geom_raster(aes(fill = reject.beta.x)) + ylab("Autoregression in Y") + xlab("Autoregression in X") + ggtitle("T=1000, Static") + guides(fill=guide_legend(title="Proportion \n Rejected")) + scale_fill_continuous(limits=c(0, .7), breaks=seq(0, .7, by = 0.1), low = "white", high = "black") + geom_contour(aes(z = reject.beta.x), breaks = prob) + geom_dl(aes(label = ..level..), method="bottom.pieces", stat="contour", breaks = prob)

pdf("yi0-xi0-SR.pdf", width = 1.618*7, height = 7.5)
grid.arrange(p1, p2, p3, p4, p5, p6, p7, p8, p9, ncol = 3)
dev.off()

# Long-run:
p1 <- ggplot(data = subset(container.collapse.ardl.50, e.var == 1), aes(x = ar.param.X, y = ar.param.Y, z = reject.lrm)) + theme_minimal() + geom_raster(aes(fill = reject.lrm)) + ylab("Autoregression in Y") + xlab("Autoregression in X") + ggtitle("T=50, ARDL/ECM") + guides(fill=guide_legend(title="Proportion \n Rejected")) + scale_fill_continuous(limits=c(0, .15), breaks=seq(0, .15, by = 0.03), low = "white", high = "black") + geom_contour(aes(z = reject.lrm), breaks = prob) + geom_dl(aes(label = ..level..), method="bottom.pieces", stat="contour", breaks = prob)

p2 <- ggplot(data = subset(container.collapse.ldv.50, e.var == 1), aes(x = ar.param.X, y = ar.param.Y, z = reject.lrm)) + theme_minimal() + geom_raster(aes(fill = reject.lrm)) + ylab("Autoregression in Y") + xlab("Autoregression in X") + ggtitle("T=50, LDV") + guides(fill=guide_legend(title="Proportion \n Rejected")) + scale_fill_continuous(limits=c(0, .15), breaks=seq(0, .15, by = 0.03), low = "white", high = "black") + geom_contour(aes(z = reject.lrm), breaks = prob) + geom_dl(aes(label = ..level..), method="bottom.pieces", stat="contour", breaks = prob)

p3 <- ggplot(data = subset(container.collapse.ardl.250, e.var == 1), aes(x = ar.param.X, y = ar.param.Y, z = reject.lrm)) + theme_minimal() + geom_raster(aes(fill = reject.lrm)) + ylab("Autoregression in Y") + xlab("Autoregression in X") + ggtitle("T=250, ARDL/ECM") + guides(fill=guide_legend(title="Proportion \n Rejected")) + scale_fill_continuous(limits=c(0, .15), breaks=seq(0, .15, by = 0.03), low = "white", high = "black") + geom_contour(aes(z = reject.lrm), breaks = prob) + geom_dl(aes(label = ..level..), method="bottom.pieces", stat="contour", breaks = prob)

p4 <- ggplot(data = subset(container.collapse.ldv.250, e.var == 1), aes(x = ar.param.X, y = ar.param.Y, z = reject.lrm)) + theme_minimal() + geom_raster(aes(fill = reject.lrm)) + ylab("Autoregression in Y") + xlab("Autoregression in X") + ggtitle("T=250, LDV") + guides(fill=guide_legend(title="Proportion \n Rejected")) + scale_fill_continuous(limits=c(0, .15), breaks=seq(0, .15, by = 0.03), low = "white", high = "black") + geom_contour(aes(z = reject.lrm), breaks = prob) + geom_dl(aes(label = ..level..), method="bottom.pieces", stat="contour", breaks = prob)

p5 <- ggplot(data = subset(container.collapse.ardl.1000, e.var == 1), aes(x = ar.param.X, y = ar.param.Y, z = reject.lrm)) + theme_minimal() + geom_raster(aes(fill = reject.lrm)) + ylab("Autoregression in Y") + xlab("Autoregression in X") + ggtitle("T=1000, ARDL/ECM") + guides(fill=guide_legend(title="Proportion \n Rejected")) + scale_fill_continuous(limits=c(0, .15), breaks=seq(0, .15, by = 0.03), low = "white", high = "black") + geom_contour(aes(z = reject.lrm), breaks = prob) + geom_dl(aes(label = ..level..), method="bottom.pieces", stat="contour", breaks = prob)

p6 <- ggplot(data = subset(container.collapse.ldv.1000, e.var == 1), aes(x = ar.param.X, y = ar.param.Y, z = reject.lrm)) + theme_minimal() + geom_raster(aes(fill = reject.lrm)) + ylab("Autoregression in Y") + xlab("Autoregression in X") + ggtitle("T=1000, LDV") + guides(fill=guide_legend(title="Proportion \n Rejected")) + scale_fill_continuous(limits=c(0, .15), breaks=seq(0, .15, by = 0.03), low = "white", high = "black") + geom_contour(aes(z = reject.beta.x), breaks = prob) + geom_dl(aes(label = ..level..), method="bottom.pieces", stat="contour", breaks = prob)

pdf("yi0-xi0-LR.pdf", width = 1.618*7, height = 7)
grid.arrange(p1, p2, p3, p4, p5, p6, ncol = 2)
dev.off()

# Short-run, increased error variance:
p1 <- ggplot(data = subset(container.collapse.ardl.50, e.var == 5), aes(x = ar.param.X, y = ar.param.Y, z = reject.beta.x)) + theme_minimal() + geom_raster(aes(fill = reject.beta.x)) + ylab("Autoregression in Y") + xlab("Autoregression in X") + ggtitle("T=50, ARDL/ECM") + guides(fill=guide_legend(title="Proportion \n Rejected")) + scale_fill_continuous(limits=c(0, .7), breaks=seq(0, .7, by = 0.1), low = "white", high = "black") + geom_contour(aes(z = reject.beta.x), breaks = prob) + geom_dl(aes(label = ..level..), method="bottom.pieces", stat="contour", breaks = prob)

p2 <- ggplot(data = subset(container.collapse.ldv.50, e.var == 5), aes(x = ar.param.X, y = ar.param.Y, z = reject.beta.x)) + theme_minimal() + geom_raster(aes(fill = reject.beta.x)) + ylab("Autoregression in Y") + xlab("Autoregression in X") + ggtitle("T=50, LDV") + guides(fill=guide_legend(title="Proportion \n Rejected")) + scale_fill_continuous(limits=c(0, .7), breaks=seq(0, .7, by = 0.1), low = "white", high = "black") + geom_contour(aes(z = reject.beta.x), breaks = prob) + geom_dl(aes(label = ..level..), method="bottom.pieces", stat="contour", breaks = prob)

p3 <- ggplot(data = subset(container.collapse.static.50, e.var == 5), aes(x = ar.param.X, y = ar.param.Y, z = reject.beta.x)) + theme_minimal() + geom_raster(aes(fill = reject.beta.x)) + ylab("Autoregression in Y") + xlab("Autoregression in X") + ggtitle("T=50, Static") + guides(fill=guide_legend(title="Proportion \n Rejected")) + scale_fill_continuous(limits=c(0, .7), breaks=seq(0, .7, by = 0.1), low = "white", high = "black") + geom_contour(aes(z = reject.beta.x), breaks = prob) + geom_dl(aes(label = ..level..), method="bottom.pieces", stat="contour", breaks = prob)

p4 <- ggplot(data = subset(container.collapse.ardl.250, e.var == 5), aes(x = ar.param.X, y = ar.param.Y, z = reject.beta.x)) + theme_minimal() + geom_raster(aes(fill = reject.beta.x)) + ylab("Autoregression in Y") + xlab("Autoregression in X") + ggtitle("T=250, ARDL/ECM") + guides(fill=guide_legend(title="Proportion \n Rejected")) + scale_fill_continuous(limits=c(0, .7), breaks=seq(0, .7, by = 0.1), low = "white", high = "black") + geom_contour(aes(z = reject.beta.x), breaks = prob) + geom_dl(aes(label = ..level..), method="bottom.pieces", stat="contour", breaks = prob)

p5 <- ggplot(data = subset(container.collapse.ldv.250, e.var == 5), aes(x = ar.param.X, y = ar.param.Y, z = reject.beta.x)) + theme_minimal() + geom_raster(aes(fill = reject.beta.x)) + ylab("Autoregression in Y") + xlab("Autoregression in X") + ggtitle("T=250, LDV") + guides(fill=guide_legend(title="Proportion \n Rejected")) + scale_fill_continuous(limits=c(0, .7), breaks=seq(0, .7, by = 0.1), low = "white", high = "black") + geom_contour(aes(z = reject.beta.x), breaks = prob) + geom_dl(aes(label = ..level..), method="bottom.pieces", stat="contour", breaks = prob)

p6 <- ggplot(data = subset(container.collapse.static.250, e.var == 5), aes(x = ar.param.X, y = ar.param.Y, z = reject.beta.x)) + theme_minimal() + geom_raster(aes(fill = reject.beta.x)) + ylab("Autoregression in Y") + xlab("Autoregression in X") + ggtitle("T=250, Static") + guides(fill=guide_legend(title="Proportion \n Rejected")) + scale_fill_continuous(limits=c(0, .7), breaks=seq(0, .7, by = 0.1), low = "white", high = "black") + geom_contour(aes(z = reject.beta.x), breaks = prob) + geom_dl(aes(label = ..level..), method="bottom.pieces", stat="contour", breaks = prob)

p7 <- ggplot(data = subset(container.collapse.ardl.1000, e.var == 5), aes(x = ar.param.X, y = ar.param.Y, z = reject.beta.x)) + theme_minimal() + geom_raster(aes(fill = reject.beta.x)) + ylab("Autoregression in Y") + xlab("Autoregression in X") + ggtitle("T=1000, ARDL/ECM") + guides(fill=guide_legend(title="Proportion \n Rejected")) + scale_fill_continuous(limits=c(0, .7), breaks=seq(0, .7, by = 0.1), low = "white", high = "black") + geom_contour(aes(z = reject.beta.x), breaks = prob) + geom_dl(aes(label = ..level..), method="bottom.pieces", stat="contour", breaks = prob)

p8 <- ggplot(data = subset(container.collapse.ldv.1000, e.var == 5), aes(x = ar.param.X, y = ar.param.Y, z = reject.beta.x)) + theme_minimal() + geom_raster(aes(fill = reject.beta.x)) + ylab("Autoregression in Y") + xlab("Autoregression in X") + ggtitle("T=1000, LDV") + guides(fill=guide_legend(title="Proportion \n Rejected")) + scale_fill_continuous(limits=c(0, .7), breaks=seq(0, .7, by = 0.1), low = "white", high = "black") + geom_contour(aes(z = reject.beta.x), breaks = prob) + geom_dl(aes(label = ..level..), method="bottom.pieces", stat="contour", breaks = prob)

p9 <- ggplot(data = subset(container.collapse.static.1000, e.var == 5), aes(x = ar.param.X, y = ar.param.Y, z = reject.beta.x)) + theme_minimal() + geom_raster(aes(fill = reject.beta.x)) + ylab("Autoregression in Y") + xlab("Autoregression in X") + ggtitle("T=1000, Static") + guides(fill=guide_legend(title="Proportion \n Rejected")) + scale_fill_continuous(limits=c(0, .7), breaks=seq(0, .7, by = 0.1), low = "white", high = "black") + geom_contour(aes(z = reject.beta.x), breaks = prob) + geom_dl(aes(label = ..level..), method="bottom.pieces", stat="contour", breaks = prob)

pdf("yi0-xi0-SR-evar5.pdf", width = 1.618*7, height = 7.5)
grid.arrange(p1, p2, p3, p4, p5, p6, p7, p8, p9, ncol = 3)
dev.off()

# Long-run, increased error variance:
p1 <- ggplot(data = subset(container.collapse.ardl.50, e.var == 5), aes(x = ar.param.X, y = ar.param.Y, z = reject.lrm)) + theme_minimal() + geom_raster(aes(fill = reject.lrm)) + ylab("Autoregression in Y") + xlab("Autoregression in X") + ggtitle("T=50, ARDL/ECM") + guides(fill=guide_legend(title="Proportion \n Rejected")) + scale_fill_continuous(limits=c(0, .15), breaks=seq(0, .15, by = 0.03), low = "white", high = "black") + geom_contour(aes(z = reject.lrm), breaks = prob) + geom_dl(aes(label = ..level..), method="bottom.pieces", stat="contour", breaks = prob)

p2 <- ggplot(data = subset(container.collapse.ldv.50, e.var == 5), aes(x = ar.param.X, y = ar.param.Y, z = reject.lrm)) + theme_minimal() + geom_raster(aes(fill = reject.lrm)) + ylab("Autoregression in Y") + xlab("Autoregression in X") + ggtitle("T=50, LDV") + guides(fill=guide_legend(title="Proportion \n Rejected")) + scale_fill_continuous(limits=c(0, .15), breaks=seq(0, .15, by = 0.03), low = "white", high = "black") + geom_contour(aes(z = reject.lrm), breaks = prob) + geom_dl(aes(label = ..level..), method="bottom.pieces", stat="contour", breaks = prob)

p3 <- ggplot(data = subset(container.collapse.ardl.250, e.var == 5), aes(x = ar.param.X, y = ar.param.Y, z = reject.lrm)) + theme_minimal() + geom_raster(aes(fill = reject.lrm)) + ylab("Autoregression in Y") + xlab("Autoregression in X") + ggtitle("T=250, ARDL/ECM") + guides(fill=guide_legend(title="Proportion \n Rejected")) + scale_fill_continuous(limits=c(0, .15), breaks=seq(0, .15, by = 0.03), low = "white", high = "black") + geom_contour(aes(z = reject.lrm), breaks = prob) + geom_dl(aes(label = ..level..), method="bottom.pieces", stat="contour", breaks = prob)

p4 <- ggplot(data = subset(container.collapse.ldv.250, e.var == 5), aes(x = ar.param.X, y = ar.param.Y, z = reject.lrm)) + theme_minimal() + geom_raster(aes(fill = reject.lrm)) + ylab("Autoregression in Y") + xlab("Autoregression in X") + ggtitle("T=250, LDV") + guides(fill=guide_legend(title="Proportion \n Rejected")) + scale_fill_continuous(limits=c(0, .15), breaks=seq(0, .15, by = 0.03), low = "white", high = "black") + geom_contour(aes(z = reject.lrm), breaks = prob) + geom_dl(aes(label = ..level..), method="bottom.pieces", stat="contour", breaks = prob)

p5 <- ggplot(data = subset(container.collapse.ardl.1000, e.var == 5), aes(x = ar.param.X, y = ar.param.Y, z = reject.lrm)) + theme_minimal() + geom_raster(aes(fill = reject.lrm)) + ylab("Autoregression in Y") + xlab("Autoregression in X") + ggtitle("T=1000, ARDL/ECM") + guides(fill=guide_legend(title="Proportion \n Rejected")) + scale_fill_continuous(limits=c(0, .15), breaks=seq(0, .15, by = 0.03), low = "white", high = "black") + geom_contour(aes(z = reject.lrm), breaks = prob) + geom_dl(aes(label = ..level..), method="bottom.pieces", stat="contour", breaks = prob)

p6 <- ggplot(data = subset(container.collapse.ldv.1000, e.var == 5), aes(x = ar.param.X, y = ar.param.Y, z = reject.lrm)) + theme_minimal() + geom_raster(aes(fill = reject.lrm)) + ylab("Autoregression in Y") + xlab("Autoregression in X") + ggtitle("T=1000, LDV") + guides(fill=guide_legend(title="Proportion \n Rejected")) + scale_fill_continuous(limits=c(0, .15), breaks=seq(0, .15, by = 0.03), low = "white", high = "black") + geom_contour(aes(z = reject.lrm), breaks = prob) + geom_dl(aes(label = ..level..), method="bottom.pieces", stat="contour", breaks = prob)

pdf("yi0-xi0-LR-evar5.pdf", width = 1.618*7, height = 7)
grid.arrange(p1, p2, p3, p4, p5, p6, ncol = 2)
dev.off()

# -- mean square error plots for SR effects, sigma_e^2 = 1
# plausible ranges for breaks:
summary(subset(container.collapse.static, e.var == 1, select = mse.beta.x))
p1 <- ggplot(data = subset(container.collapse.ardl.50, e.var == 1), aes(x = ar.param.X, y = ar.param.Y, z = mse.beta.x)) + theme_minimal() + geom_raster(aes(fill = mse.beta.x)) + ylab("Autoregression in Y") + xlab("Autoregression in X") + ggtitle("T=50, ARDL/ECM") + guides(fill=guide_legend(title="MSE")) + scale_fill_continuous(limits=c(0, .3), breaks=seq(0, .3, by = 0.05), low = "white", high = "black")

p2 <- ggplot(data = subset(container.collapse.ldv.50, e.var == 1), aes(x = ar.param.X, y = ar.param.Y, z = mse.beta.x)) + theme_minimal() + geom_raster(aes(fill = mse.beta.x)) + ylab("Autoregression in Y") + xlab("Autoregression in X") + ggtitle("T=50, LDV") + guides(fill=guide_legend(title="MSE")) + scale_fill_continuous(limits=c(0, .3), breaks=seq(0, .3, by = 0.05), low = "white", high = "black")

p3 <- ggplot(data = subset(container.collapse.static.50, e.var == 1), aes(x = ar.param.X, y = ar.param.Y, z = mse.beta.x)) + theme_minimal() + geom_raster(aes(fill = mse.beta.x)) + ylab("Autoregression in Y") + xlab("Autoregression in X") + ggtitle("T=50, Static") + guides(fill=guide_legend(title="MSE")) + scale_fill_continuous(limits=c(0, .3), breaks=seq(0, .3, by = 0.05), low = "white", high = "black")

p4 <- ggplot(data = subset(container.collapse.ardl.250, e.var == 1), aes(x = ar.param.X, y = ar.param.Y, z = mse.beta.x)) + theme_minimal() + geom_raster(aes(fill = mse.beta.x)) + ylab("Autoregression in Y") + xlab("Autoregression in X") + ggtitle("T=250, ARDL/ECM") + guides(fill=guide_legend(title="MSE")) + scale_fill_continuous(limits=c(0, .3), breaks=seq(0, .3, by = 0.05), low = "white", high = "black")

p5 <- ggplot(data = subset(container.collapse.ldv.250, e.var == 1), aes(x = ar.param.X, y = ar.param.Y, z = mse.beta.x)) + theme_minimal() + geom_raster(aes(fill = mse.beta.x)) + ylab("Autoregression in Y") + xlab("Autoregression in X") + ggtitle("T=250, LDV") + guides(fill=guide_legend(title="MSE")) + scale_fill_continuous(limits=c(0, .3), breaks=seq(0, .3, by = 0.05), low = "white", high = "black")

p6 <- ggplot(data = subset(container.collapse.static.250, e.var == 1), aes(x = ar.param.X, y = ar.param.Y, z = mse.beta.x)) + theme_minimal() + geom_raster(aes(fill = mse.beta.x)) + ylab("Autoregression in Y") + xlab("Autoregression in X") + ggtitle("T=250, Static") + guides(fill=guide_legend(title="MSE")) + scale_fill_continuous(limits=c(0, .3), breaks=seq(0, .3, by = 0.05), low = "white", high = "black")

p7 <- ggplot(data = subset(container.collapse.ardl.1000, e.var == 1), aes(x = ar.param.X, y = ar.param.Y, z = mse.beta.x)) + theme_minimal() + geom_raster(aes(fill = mse.beta.x)) + ylab("Autoregression in Y") + xlab("Autoregression in X") + ggtitle("T=1000, ARDL/ECM") + guides(fill=guide_legend(title="MSE")) + scale_fill_continuous(limits=c(0, .3), breaks=seq(0, .3, by = 0.05), low = "white", high = "black")

p8 <- ggplot(data = subset(container.collapse.ldv.1000, e.var == 1), aes(x = ar.param.X, y = ar.param.Y, z = mse.beta.x)) + theme_minimal() + geom_raster(aes(fill = mse.beta.x)) + ylab("Autoregression in Y") + xlab("Autoregression in X") + ggtitle("T=1000, LDV") + guides(fill=guide_legend(title="MSE")) + scale_fill_continuous(limits=c(0, .3), breaks=seq(0, .3, by = 0.05), low = "white", high = "black")

p9 <- ggplot(data = subset(container.collapse.static.1000, e.var == 1), aes(x = ar.param.X, y = ar.param.Y, z = mse.beta.x)) + theme_minimal() + geom_raster(aes(fill = mse.beta.x)) + ylab("Autoregression in Y") + xlab("Autoregression in X") + ggtitle("T=1000, Static") + guides(fill=guide_legend(title="MSE")) + scale_fill_continuous(limits=c(0, .3), breaks=seq(0, .3, by = 0.05), low = "white", high = "black")

pdf("yi0-xi0-SR-MSE.pdf", width = 1.618*7, height = 7.5)
grid.arrange(p1, p2, p3, p4, p5, p6, p7, p8, p9, ncol = 3)
dev.off()

# -- median square error plots for LR effects, sigma_e^2 = 1
summary(subset(container.collapse.ardl, e.var == 1, select = mse.lrm))
summary(subset(container.collapse.ecm, e.var == 1, select = mse.lrm))
summary(subset(container.collapse.ldv, e.var == 1, select = mse.lrm))

p1 <- ggplot(data = subset(container.collapse.ardl.50, e.var == 1), aes(x = ar.param.X, y = ar.param.Y, z = mse.lrm)) + theme_minimal() + geom_raster(aes(fill = mse.lrm)) + ylab("Autoregression in Y") + xlab("Autoregression in X") + ggtitle("T=50, ARDL/ECM") + guides(fill=guide_legend(title="Med Sq Err")) + scale_fill_continuous(limits=c(0, 1.5), breaks=seq(0, 1.5, by = .3), low = "white", high = "black")

p2 <- ggplot(data = subset(container.collapse.ldv.50, e.var == 1), aes(x = ar.param.X, y = ar.param.Y, z = mse.lrm)) + theme_minimal() + geom_raster(aes(fill = mse.lrm)) + ylab("Autoregression in Y") + xlab("Autoregression in X") + ggtitle("T=50, LDV") + guides(fill=guide_legend(title="Med Sq Err")) + scale_fill_continuous(limits=c(0, 1.5), breaks=seq(0, 1.5, by = .3), low = "white", high = "black")

p3 <- ggplot(data = subset(container.collapse.ardl.250, e.var == 1), aes(x = ar.param.X, y = ar.param.Y, z = mse.lrm)) + theme_minimal() + geom_raster(aes(fill = mse.lrm)) + ylab("Autoregression in Y") + xlab("Autoregression in X") + ggtitle("T=250, ARDL/ECM") + guides(fill=guide_legend(title="Med Sq Err")) + scale_fill_continuous(limits=c(0, 1.5), breaks=seq(0, 1.5, by = .3), low = "white", high = "black")

p4 <- ggplot(data = subset(container.collapse.ldv.250, e.var == 1), aes(x = ar.param.X, y = ar.param.Y, z = mse.lrm)) + theme_minimal() + geom_raster(aes(fill = mse.lrm)) + ylab("Autoregression in Y") + xlab("Autoregression in X") + ggtitle("T=250, LDV") + guides(fill=guide_legend(title="Med Sq Err")) + scale_fill_continuous(limits=c(0, 1.5), breaks=seq(0, 1.5, by = .3), low = "white", high = "black")

p5 <- ggplot(data = subset(container.collapse.ardl.1000, e.var == 1), aes(x = ar.param.X, y = ar.param.Y, z = mse.lrm)) + theme_minimal() + geom_raster(aes(fill = mse.lrm)) + ylab("Autoregression in Y") + xlab("Autoregression in X") + ggtitle("T=1000, ARDL/ECM") + guides(fill=guide_legend(title="Med Sq Err")) + scale_fill_continuous(limits=c(0, 1.5), breaks=seq(0, 1.5, by = .3), low = "white", high = "black")

p6 <- ggplot(data = subset(container.collapse.ldv.1000, e.var == 1), aes(x = ar.param.X, y = ar.param.Y, z = mse.lrm)) + theme_minimal() + geom_raster(aes(fill = mse.lrm)) + ylab("Autoregression in Y") + xlab("Autoregression in X") + ggtitle("T=1000, LDV") + guides(fill=guide_legend(title="Med Sq Err")) + scale_fill_continuous(limits=c(0, 1.5), breaks=seq(0, 1.5, by = .3), low = "white", high = "black")

pdf("yi0-xi0-LR-MSE.pdf", width = 1.618*7, height = 7)
grid.arrange(p1, p2, p3, p4, p5, p6, ncol = 2)
dev.off()

# -- mean square error plots for SR effects, sigma_e^2 = 5
# plausible ranges for breaks:
summary(subset(container.collapse.static, e.var == 5, select = mse.beta.x))
p1 <- ggplot(data = subset(container.collapse.ardl.50, e.var == 5), aes(x = ar.param.X, y = ar.param.Y, z = mse.beta.x)) + theme_minimal() + geom_raster(aes(fill = mse.beta.x)) + ylab("Autoregression in Y") + xlab("Autoregression in X") + ggtitle("T=50, ARDL/ECM") + guides(fill=guide_legend(title="MSE")) + scale_fill_continuous(limits=c(0, 1.4), breaks=seq(0, 1.4, by = .2), low = "white", high = "black")

p2 <- ggplot(data = subset(container.collapse.ldv.50, e.var == 5), aes(x = ar.param.X, y = ar.param.Y, z = mse.beta.x)) + theme_minimal() + geom_raster(aes(fill = mse.beta.x)) + ylab("Autoregression in Y") + xlab("Autoregression in X") + ggtitle("T=50, LDV") + guides(fill=guide_legend(title="MSE")) + scale_fill_continuous(limits=c(0, 1.4), breaks=seq(0, 1.4, by = .2), low = "white", high = "black")

p3 <- ggplot(data = subset(container.collapse.static.50, e.var == 5), aes(x = ar.param.X, y = ar.param.Y, z = mse.beta.x)) + theme_minimal() + geom_raster(aes(fill = mse.beta.x)) + ylab("Autoregression in Y") + xlab("Autoregression in X") + ggtitle("T=50, Static") + guides(fill=guide_legend(title="MSE")) + scale_fill_continuous(limits=c(0, 1.4), breaks=seq(0, 1.4, by = .2), low = "white", high = "black")

p4 <- ggplot(data = subset(container.collapse.ardl.250, e.var == 5), aes(x = ar.param.X, y = ar.param.Y, z = mse.beta.x)) + theme_minimal() + geom_raster(aes(fill = mse.beta.x)) + ylab("Autoregression in Y") + xlab("Autoregression in X") + ggtitle("T=250, ARDL/ECM") + guides(fill=guide_legend(title="MSE")) + scale_fill_continuous(limits=c(0, 1.4), breaks=seq(0, 1.4, by = .2), low = "white", high = "black")

p5 <- ggplot(data = subset(container.collapse.ldv.250, e.var == 5), aes(x = ar.param.X, y = ar.param.Y, z = mse.beta.x)) + theme_minimal() + geom_raster(aes(fill = mse.beta.x)) + ylab("Autoregression in Y") + xlab("Autoregression in X") + ggtitle("T=250, LDV") + guides(fill=guide_legend(title="MSE")) + scale_fill_continuous(limits=c(0, 1.4), breaks=seq(0, 1.4, by = .2), low = "white", high = "black")

p6 <- ggplot(data = subset(container.collapse.static.250, e.var == 5), aes(x = ar.param.X, y = ar.param.Y, z = mse.beta.x)) + theme_minimal() + geom_raster(aes(fill = mse.beta.x)) + ylab("Autoregression in Y") + xlab("Autoregression in X") + ggtitle("T=250, Static") + guides(fill=guide_legend(title="MSE")) + scale_fill_continuous(limits=c(0, 1.4), breaks=seq(0, 1.4, by = .2), low = "white", high = "black")

p7 <- ggplot(data = subset(container.collapse.ardl.1000, e.var == 5), aes(x = ar.param.X, y = ar.param.Y, z = mse.beta.x)) + theme_minimal() + geom_raster(aes(fill = mse.beta.x)) + ylab("Autoregression in Y") + xlab("Autoregression in X") + ggtitle("T=1000, ARDL/ECM") + guides(fill=guide_legend(title="MSE")) + scale_fill_continuous(limits=c(0, 1.4), breaks=seq(0, 1.4, by = .2), low = "white", high = "black")

p8 <- ggplot(data = subset(container.collapse.ldv.1000, e.var == 5), aes(x = ar.param.X, y = ar.param.Y, z = mse.beta.x)) + theme_minimal() + geom_raster(aes(fill = mse.beta.x)) + ylab("Autoregression in Y") + xlab("Autoregression in X") + ggtitle("T=1000, LDV") + guides(fill=guide_legend(title="MSE")) + scale_fill_continuous(limits=c(0, 1.4), breaks=seq(0, 1.4, by = .2), low = "white", high = "black")

p9 <- ggplot(data = subset(container.collapse.static.1000, e.var == 5), aes(x = ar.param.X, y = ar.param.Y, z = mse.beta.x)) + theme_minimal() + geom_raster(aes(fill = mse.beta.x)) + ylab("Autoregression in Y") + xlab("Autoregression in X") + ggtitle("T=1000, Static") + guides(fill=guide_legend(title="MSE")) + scale_fill_continuous(limits=c(0, 1.4), breaks=seq(0, 1.4, by = .2), low = "white", high = "black")

pdf("yi0-xi0-SR-MSE-evar5.pdf", width = 1.618*7, height = 7.5)
grid.arrange(p1, p2, p3, p4, p5, p6, p7, p8, p9, ncol = 3)
dev.off()


# -- median square error plots for LR effects, sigma_e^2 = 5
summary(subset(container.collapse.ardl, e.var == 5, select = mse.lrm))
p1 <- ggplot(data = subset(container.collapse.ardl.50, e.var == 5), aes(x = ar.param.X, y = ar.param.Y, z = mse.lrm)) + theme_minimal() + geom_raster(aes(fill = mse.lrm)) + ylab("Autoregression in Y") + xlab("Autoregression in X") + ggtitle("T=50, ARDL/ECM") + guides(fill=guide_legend(title="Med Sq Err")) + scale_fill_continuous(limits=c(0, 7), breaks=seq(0, 7, by = 1), low = "white", high = "black")

p2 <- ggplot(data = subset(container.collapse.ldv.50, e.var == 5), aes(x = ar.param.X, y = ar.param.Y, z = mse.lrm)) + theme_minimal() + geom_raster(aes(fill = mse.lrm)) + ylab("Autoregression in Y") + xlab("Autoregression in X") + ggtitle("T=50, LDV") + guides(fill=guide_legend(title="Med Sq Err")) + scale_fill_continuous(limits=c(0, 7), breaks=seq(0, 7, by = 1), low = "white", high = "black")

p3 <- ggplot(data = subset(container.collapse.ardl.250, e.var == 5), aes(x = ar.param.X, y = ar.param.Y, z = mse.lrm)) + theme_minimal() + geom_raster(aes(fill = mse.lrm)) + ylab("Autoregression in Y") + xlab("Autoregression in X") + ggtitle("T=250, ARDL/ECM") + guides(fill=guide_legend(title="Med Sq Err")) + scale_fill_continuous(limits=c(0, 7), breaks=seq(0, 7, by = 1), low = "white", high = "black")

p4 <- ggplot(data = subset(container.collapse.ldv.250, e.var == 5), aes(x = ar.param.X, y = ar.param.Y, z = mse.lrm)) + theme_minimal() + geom_raster(aes(fill = mse.lrm)) + ylab("Autoregression in Y") + xlab("Autoregression in X") + ggtitle("T=250, LDV") + guides(fill=guide_legend(title="Med Sq Err")) + scale_fill_continuous(limits=c(0, 7), breaks=seq(0, 7, by = 1), low = "white", high = "black")

p5 <- ggplot(data = subset(container.collapse.ardl.1000, e.var == 5), aes(x = ar.param.X, y = ar.param.Y, z = mse.lrm)) + theme_minimal() + geom_raster(aes(fill = mse.lrm)) + ylab("Autoregression in Y") + xlab("Autoregression in X") + ggtitle("T=1000, ARDL/ECM") + guides(fill=guide_legend(title="Med Sq Err")) + scale_fill_continuous(limits=c(0, 7), breaks=seq(0, 7, by = 1), low = "white", high = "black")

p6 <- ggplot(data = subset(container.collapse.ldv.1000, e.var == 5), aes(x = ar.param.X, y = ar.param.Y, z = mse.lrm)) + theme_minimal() + geom_raster(aes(fill = mse.lrm)) + ylab("Autoregression in Y") + xlab("Autoregression in X") + ggtitle("T=1000, LDV") + guides(fill=guide_legend(title="Med Sq Err")) + scale_fill_continuous(limits=c(0, 7), breaks=seq(0, 7, by = 1), low = "white", high = "black")

pdf("yi0-xi0-LR-MSE-evar5.pdf", width = 1.618*7, height = 7)
grid.arrange(p1, p2, p3, p4, p5, p6, ncol = 2)
dev.off()